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ABSTRACT 

In this paper we consider the effect of the direct ionizing stellar radiation field on 
the evolution of protoplanetary discs subject to photoevaporative winds. We suggest 
that models which com bine viscous evolution with photoevaporation of the disc (e.g. 
IClarke. Gendrin & Sotomavor 2001) incorrectly neglect the direct field after the inner disc 
has drained, at late times in the evolution. We construct models of the photoevaporative wind 
produced by the direct field, first using simple analytic arguments and later using detailed 
numerical hydrodynamics. We find that the wind produced by the direct field at late times is 
much larger than has previously been assumed, and we show that the mass-loss rate scales as 

1 /2 

(where is the radius of the instantaneous inner disc edge). We suggest that this result 
has important consequences for theories of disc evolution, and go on to consider the effects 
of this result on disc evolution in detail in a companion paper (lAlexander. Clarke & Pringl3 
|2006b). 

Key words: accretion, accretion discs - circumstellar matter - hydrodynamics - planetary 
systems: protoplanetary discs - stars: pre-main-sequence 



; 1 INTRODUCTION 

The evolution and eventual dispersal of discs around young stars 
is an important area of study, for theories of both star and planet 
formation. It is now well-established that at an age of ~ lO^yr 
most stars are surrounded by discs t hat are optically thick at op- 
tical and infrared wavelengths (e.g. iKenvon & Hartmarmlll99^ 
and typically have masses of a few percent of a solar mass 
(e.g. Beckwith et al. 1990). However at an age of ~ lO^yr most 
stars are not seen to have di scs, suggesting that disc lifetimes 
are typically a few Myr (e.g. lHaisch. Lada & Ladai.2001.) . How 
stars lose their discs remains an unsolved question. Observations 
of T Tauri stars (TTs) in the infrared and at millimetre wave- 
lengths detect few objects intermediate between the disc-bearing 
classi cal T Tauri (CTT) and disc-less weak-lined T Tauri (WTT) 
states ISkrutskie et al.*199dl:lKenvo n & Hartmannll99^|Persi et alJ 

[2000 ; Andrews & Williams 2005). This suggests that the transition 
from CTT to WTT is very rapid, with the dispersal time estimated 
to be ~ lOVr 1 Simon & Prato 1995; Wolk & Walter 1996). More- 
over the simultaneous decline is disc emission across a wide range 
in wavelength suggests that the dispersal is essentially simultane- 
ous across the entire radial extent of the disc. 

This "two-time-scale" behaviour is inconsistent with con- 
ventional models of disc evolution, which predict power- 
law declines in the disc and theref ore predict dispersal 
times comparable to the disc lifetimes iHartmann etalJ 1199a 

lArmitage. Clarke & Tout 1999.) . However it has been shown that 



models which combine photoevaporation of the disc material with 
viscous evolution can reproduce th is two-time-scale behaviour 
iClarke. Gendrin & SotomavojEoOll) . Models of so-called pho- 
toevapor ative disc winds were first suggested at least 20 years 
ago te.g. lBallv & Scovilldll98j) . However the first models which 
treated the flow and radiative transfer in a fully self-consistent 
manner were repo rted by Hollenbach, lohnstone & Shu ( 1993) and 
LShu, lohnstone & Hollenbach ( 1993), and then extended in depth 
by ^ollenbach e t aL U994) . Further studies have extended the 
study o f the detai ls in a number of ways (e.g. the effects of dust, 
iRichling & York3ll997t the inclusion of external radiation fields 
and n on-ionizing UV radiation, Johns t Qne,JIgllenb ach & Balfvl 
[l998; detailed numerical hydrodynamics . Ipont et al 120041) . but the 
underlying principles remain the same. Ionizing radiation from the 
central star produces a hot ionized layer on the surface of the disc, 
with conditions akin to an Hii region. Near to the star the ion- 
ized gas remains bound to the disc, but outside some critical radius 
(known as the gravitational radius and denoted by Rg) the ionized 
layer becomes unbound and flows from the disc surface as a wind. 
For solar-mass stars such as TTs the gravitational radius is typi- 
cally a few AU. For a constant ionizing flux O the wind rate is con- 
stant, as long as the disc remains optically thick to ionizing photons 
perpendicular to the plane of the disc. Furthermore the diffuse (re- 
combination) radiation field is the dominant source of ionizations 
at all radii of interest'. In their "weak-wind" case, applicable to 
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1 iHoUenbach eTS] 1 1994 found that the radius beyond which the difl'use 
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TTs. lHollenbach et alj il99il found that the mass-loss rate due to 
photoevaporation is given by 



10 km/s flow 



: 4.4 X 10" 



10"' s- 



Moyr 



(1) 



More recently, hydrodynamic modelling of a photoevaporative 
wind has resulted in slight modification of this result. When hy- 
drodynami c effect are considered the "effec tive /?g" is reduced by a 
factor of 5 iLiffmarJ200.ilFont et alllOol. and the mass-loss rate 
is reduced by a factor of around 3 <Font et alj2004. However the 
qualita tive behaviour is unchanged from that of iHollenbach et alj 

The so-called "UV-switch" model of IClarke et all <200lh 
couples a photoevaporative wind to a simple disc evolution model. 
At early times the accretion rate through the disc is much larger 
than the wind rate, and the wind has a negligible effect. However 
at late times photoevaporation becomes important, depriving the 
disc of resupply inside R^. At this point the inner disc drains on 
its own, short, viscous time-scale, giving a dispersal time much 
shorter than the disc lifetime . A number of similar studies have 
now b een conducted jMatsuvama, Johnstone & HartmannI 
| 2003l lArmitase. Clarke & Fallal 120031: Rudea l2004t 
iTakeuchi. Clarke & LirJl2005h . and this class of models shows a 
number of attractive properties. 

However Clarke et al. ( 20(^ highlighted two key problems 
with the model. Firstly, the model requires that TTs produce a 
rather large ionizing flux (of order 10"" ionizing photons per sec- 
ond), and secondly that the outer disc, beyond R„ is dispersed much 
too slowly to satisfy millimetre observations of TT discs. We have 
previously shown that it is reasonable to treat TT chromospheres 
as having a constant ionizing flux in the range ~ 10'"-10'''*s"' 
jAlexander. Clarke & Pringldl2005h . We now seek to address the 
"outer disc problem" by highlighting a flaw in the original UV- 
switch model. The photoevaporative wind models described above 
are based on the premise that the disc is optically thick to Lyman 
continuum photons at all radii. Consequently the diffuse (recombi- 
nation) radiation field is the dominant source of ionizing photons 
at large radii (HoUenbach et al. 1994), and drives the photoevap- 
orative wind. However once the inner disc is drained it becomes 
optically thin to Lyman continuum photons, eliminating this diffuse 
field and rendering the simple wind prescription used above invalid. 
Instead, once the inner disc has drained we must consider direct 
ionization of the inner edge of the disc. Clarke et al. ( 2001) neglect 
this process, and find that the time-scale for dispersal of the outer 
disc is limited by the time material takes to diffuse inward to (as 
the most of the mass-loss occurs close to R ^). Further, the /?~^^^ de- 
pendence of the wind profile at large radii fHollenbach et al .'1994') 
means that the mass-loss rate due to photoevaporation decreases 
significantly with time as the inner edge of the disc moves outward. 
Consequently the dispersal of the outer disc occurs on the viscous 
time-scale of the outer disc, and thus the outer disc is dispersed in 
a time comparable to the disc lifetime, much too slowly to satisfy 
observational constraints. We suggest that photoevaporation by the 
direct radiation field r esults in a dispersal time significantly shorter 
than that predicted bv lClarke et alj i200ll) . and now seek to model 
the effects of this process on the evolution of the outer disc. 

In this paper we model the photoevaporative wind produced 
by the direct field. In Section |5| we present an analytic treatment 



field dominate.s is 3.8 (M./IMq)"^ (R./1Rq)^'' Ro, which is much smaller 
than R„ for parameters typical of TTs. 
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Figure 1. Schematic representation of the disc wind produced by direct 
photoevaporation. Ionizing radiation from the star creates an ionized layer 
on the disc surface. The ionized gas is unbound, and flows away from the 
disc at approximately the sound speed. 



of the problem, highlighting the important physical processes. In 
Section|3|we extend our analysis by using numerical hydrodynam- 
ics, and describe the numerical simulations we have performed. In 
Section|4]we present our results, and show that the wind profile is 
well-matched by a simple analytic form. In Section|5|we discuss the 
caveats that apply to our modelling, as well as the consequences of 
our results for disc evolution models, and in Section|6|we summa- 
rize our conclusions. We consider the effects of o ur results on disc 
evolution models in detail in a companion paper i Alexander et alJ 
■2006h. hereafter Paper II). 



2 ANALYTIC MODEL 

We now consider the problem which we refer to as "direct photoe- 
vaporation". As discussed above, after the inner-disc draining of the 
UV-switch model the inner disc becomes optically thin to Lyman 
continuum photons, and so we must consider the influence of the 
direct radiation field rather than the diffuse field considered previ- 
ously. We note that while the UV-switch model shows that the gas 
in the inner disc will become optically thin to ionizing radiation, 
some dust opacity may remain after the time at which the inner gas 
disc drains. An detailed treatment of gas-grain dynamics is beyond 
the scope of this work, but we will consider this issue qualitatively 
in the discussion (Section|5}. We expect similar physics to apply in 
the case of the direct field as in the diffuse field case, with the ioniz- 
ing radiation creating a thin ionized layer on the surface of the disc 
(see Fig0. However the disc is truncated at some inner radius such 
that the ionized layer is unbound, and flows away from the surface 
as a disc wind. Additionally, at the inner edge any flow perpendicu- 
lar to the disc surface will move inward and so some material may 
be accreted on to the star (although it must presumably lose angu- 
lar momentum to do so). Obviously this is a complicated dynamic 
process, and in the next section we use numerical hydrodynamics to 
model the dynamics of direct photoevaporation. However we first 
seek a theoretical framework in which to place these hydrodynamic 
models. 

In order to study the problem of direct photoevapora- 
tion analytically we us e a similar approach to that adopted by 
IHollenbach e7S] jl994b . We assume that the mass-loss per unit 
area from the disc at a given point is given by pc^, where p is the 
density at the base of the region layer and Cs is the sound speed of 
the ionized gas. Consequently the density at the at the base of the 
ionized region is critical to the determination of the photoevapora- 
tive wind. 
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Figure 2. Diagram showing geometry of tlie direct ionization problem. 
Along a line-of-sight with polar angle 6 and angular size A6, the angle 
between the ionization front and the line-of-sight is /?. Consequently the 
projected length of the column perpendicular to the front is rA0/ sin/?. 



We solve for the number density at the base of the ionized 
region (i.e. at the ionization front) as follows. We neglect recom- 
binations between the source and the disc surface, and assume that 
the location of ionization front is determined by ionization balance 
only. We also assume azimuthal symmetry, and integrate over the 
azimuthal coordinate throughout. Consequently, along any given 
line-of-sight from the source the rate of recombinations, A'rec, in a 
volume AV at the ionization front must balance the rate of ionizing 
photons absorbed at the front, N\a^. A column with polar angle Q 
and angular size A(? has a total area equal to 2nr^ sin SAS, so for an 
ionizing flux <b the ionization rate at the front is- 

A'.on = ^ sin 6)AeO . (2) 

The geometry of the problem is shown in Fig|5| We have assumed 
that at the ionization front we have ionization balance within some 
volume AV. Perturbations to the disc structure on a length-scale 
shorter than the disc scale-height H{R) are unlikely to be dynami- 
cally stable (e.g. lLin. Paoa loizou & Faulkner" 1985'\ so we assume 
that Al/ has a thickness (perpendicular to the front) that is equal to 
the disc scale-height H(K). We verify this assumption a posteriori 
in Section l43l Consequently AV is given by 



AV = InRH 



rAe 
sin/3 ' 



(3) 



where /3 is the angle between the ray-path and the ionization front 
(see Fig|2j. The spherical radius r is related to the cylindrical radius 
Rhy R = rsin 6, and so the recombination rate is 



Wi-ec = InanlRH- 



RAe 



(4) 



sin0sin/3 

where no(R) is the number density at the ionization front and a is 



At various points in this paper it is convenient to work in either cylin- 
drical or spherical polar coordinates. To avoid confusion we adopt the no- 
tation (r, d, <f>) for spherical coordinates and (S, z, ip) for cylindrical. Thus 
upper-case R always represents cylindrical radius, while lower-case r de- 
notes spherical radius. 



the hydrogen recombination coefficient. If we equate the ionization 
and recombination rates and re-arrange we find that 



/Osin-(9sin/3\''^ 



(5) 



In general, therefore, the density at the ionization front at a given 
radius R depends on the geometry of the ionization front. This in 
turn depends of the density structure of the unperturbed disc, and 
the form of the density profile does not permit exact analytic so- 
lutions. However we can look at the behaviour of the density in 
specific cases. For example, along the midplane (6 = njl) we ex- 
pect the ionization front to be perpendicular to the line-of-sight to 
the source, and so sin /3 = 1. In most discs H{R) is an increas- 
ing function of R, and at radii beyond the inner disc edge we have 
sinS < 1 and sin/? < 1, so iioiR) falls off more steeply than i?"''^ 
The exact profile of the base density is determined by the geom- 
etry of the ionization front, which must be evaluated numerically. 
However before we consider numerical modelling of the wind it is 
instructive to frame the problem in terms of known parameters and 
unknown scaling constants. 

We note at this point that the form we have derived for no(^) 
(Equation|5} is qualitatively s imilar to the "strong wind" solution 
derived by Hollenbach et al] ^993)- In this case a strong stellar 
wind modifies the disc structure out to a radius beyond Rg, "driv- 
ing down" the disc atmosphere and enabling radiation to penetrate 
to lower heights above the disc midplane at larger radii. Conse- 
quently the ;io °^ R^^^~ scaling, which applies only to the inner disc 
in the weak-wind case, extends well beyond R, in the presence of a 
strong stellar wind. Such a wind is not found for TTs, and instead 
this case is more readily applicable to O- and B-type stars. How- 
ever the geometry of the radiative transfer problem is similar to that 
considered here, with a similar lack of a disc atmosphere at small 
radii. The numerical scaling will differ, as even in the strong wind 
case the diffuse field still dominates, but we expect the form of the 
mass loss due to the direct field to be qualitatively similar to that 
derived by Hollenba ch et alJ il994) for the case of a strong stellar 
wind. 



2.1 Scaling parameters 

Following Holl enbach et alj il994h . we evaluate the wind mass- 
loss rate per unit area as 



2„i„d(^, t) = Ifim^noiR, t)ui(R, t) , 



(6) 



where /j is the mean molecular weight of the gas, mn is the mass 
of a hydrogen atom, and the factor of 2 accounts for mass-loss 
from both sides of the disc. The launch velocity U](R,t) is of or- 
der the sound speed of the ionized gas. For the moment we make 
the simplifying assumption that at the inner disc edge AV has a 
radial thickness comparable to the vertical scale-height. As noted 
above, along the midplane we have sinS = sin/? = 1, and so the 
inner edge density ni„ can be expressed as 



C 



4na(H/Rk„Rl 



1/2 



(7) 



where R^^ is the radius of the inner disc edge and C is an order- 
of-unity scaling constant. We further assume that the base density 
no(R) (Equation|5} can be related to the inner edge density «i„ by a 
dimensionless shape function f (R/R,„), which depends on the disc 
structure (and therefore depends on the form of H(R)). Therefore 
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no(R) = "in/I — 



and the launch velocity can be written as 
u,(R) = Dc, . 



(8) 



(9) 



Here D is another order-of-unity scaling constant and c, = 10km 
s"' is the sound speed of the ionized gas. Therefore the mass-loss 
profile takes the form 



2wind(^,0 = 2CDiimnCsn,Mf 



R 



R,M 



(10) 



The total mass-loss rate is found by integrating this from Ri„ to 
some outer radius i?oiit: 



M{<Roud= [ 



2nRI.„M{R)dR . 



(11) 



Substituting for ni„ from Equation0and integrating gives 



M(< Rout) = 4nCDfiiniiC\ 



1/2 



R 



1/2 



4na(H/R}„ 

r-«o.„/Ri„ 
X I xf(x)dx. 



(12) 



Typically, for TT discs, R„ut » Rm- Further, if we assume that the 
integral converges as Rout/Rm — > °° (i-e- that the form of f(x) falls 
off faster than x^-^), then we see that the mass-loss rate from the 
wind scales as RII^. Consequently we expect the mass-loss rate to 
increase as the inner edge of the disc evolves outwards. If we re- 
scale in terms of parameters typical of TTs (taking {H/R)i„ = 0.05), 
we find 



M(< Rout) = 1.74 X 10"' CZ)yU 



1/2 



\3AU/ 



1/2 



10«s-' 

•«oul/«in 

x/(x)dxMoyr"' , 



(13) 



As noted earlier, the shape function /(x) depends on both the radial 
density profile and the form of H(R). We expect a profile that falls 
off faster than R^^'^, but require numerical solution in order to de- 
termine the exact form of f(x). Moreover we must also determine 
the scaling constants C and D numerically, so we turn to numerical 
simulations in order to pursue this problem further. 



3 HYDRODYNAMIC MODEL 

In this section we construct and evaluate hydrodynamic models of 
the disc wind produced by direct photoevaporation. We first de- 
scribe the computational method used to investigate the problem 
and the simulations used to study the behaviour of the wind. We 
then present our results, which are framed in terms of the scaling 
parameters described above, before discussing the various caveats 
which apply to the models. 



3.1 Computational method 

To investigate the dynamics of direct photoevaporation we use 
the two-dimensional (2-D) grid-based hydrodynamics code zeus2d 
(Stone & Norman 1992). zeus2d uses operator splitting to integrate 
the hydrodynamic equations forward in time. Each equation is split 
into parts, which are then evaluated consecutively using the results 
from the previous evaluation. Further splitting is done between 
the two grid directions, with the order alternating on consecutive 



timesteps. The hydrodynamic equations are split into source and 
transport terms. The source terms are evaluated first, and the new 
values are then transported by computing fluxes across cell faces. 
The form of this computation makes it relatively straightforward to 
add heating due to ionizing radiation, as it is easy to update the lo- 
cal energy at the end of each timestep, before subsequent evaluation 
of the source terms. 

zeus2d provides several alternative options for the interpola- 
tion scheme and artificial viscosity. We adopt the standard van Leer 
(second-order) interpolation scheme, and the standard von Neu- 
mann & Richtmyer form for the artificial viscosity (with ^visc = 
2.0). In evaluating the timestep we adopt a Courant number of 0.4. 
Additionally, we adopt the ideal gas equation of state 



p = (r- ^)e 



(14) 



throughout, where p is gas pressure, e is energy density (per unit 
volume)', and y = 5/3 is the adiabatic exponent. 



3.1.1 Ionizing radiation 

In order to model the effects of ionizing radiation in zeus2d it is nec- 
essary to simplify the radiative transfer problem, as it is not prac- 
tical to make a full radiative transfer calculation at each timestep. 
In order to solve the problem we assume that the only significant 
heat source is the absorption of Lyman continuum photons, and 
that the only significant coolant is the re-emission of such photons 
by radiative recombination. Consequently the problem of evaluat- 
ing the heating and cooling due to radiation is reduced to a prob- 
lem of computing ionization balance. We assume that the gas is 
either extremely optically thick to ionizing photons produced by 
recombination, in which case they are absorbed locally (Case B re- 
combination), or assume that the gas is extremely optically thin to 
recombination photons, in which case they escape from the system 
(Case A). Consequently recombination photons can be neglected in 
the computation of ionization balance. This approximation (specif- 
ically the assumption of local absorption) is known as the "on-the- 
spot" (henceforth OTS) approximation. Where the OTS approxi- 
mation is not valid is the case where ionizing photons produced 
by recombinations travel some distance before being re-absorbed. 
We discuss the validity of the OTS approximation a posteriori in 
Section|5| 

In order to model the effects of a central ionizing source we 
have added two routines to the code. These routines require that a 
polar [i.e. (r,d)] coordinate grid is used, so that the radial columns 
represent ray paths from the central source. The first routine is per- 
formed at the end of each timestep (i.e. after the transport step), and 
solves for ionization balance along each ray path in order to find the 
location of the ionization front. It then sets a flag array to indicate 
whether cells should be ionized, neutral or "boundary". To do this 
we solve the equation of ionization balance along each radial grid 
column (denoted by subscript j) 



dNj 1 v-i 9 

— i = -sin0,Ae,cl,-2.^A;«^,K,, 



(15) 



where Nj is the total number of ionized atoms in the along the jth 



Note that for computational reasons zeus2d works with the energy den- 
sity per unit volume, denoted by e, rather than the more conventional energy 
per unit mass. Here e is equal to the quantity denoted by "pe" in most hy- 
drodynamic texts. 
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radial grid column, a the recombination rate coefficient, n, , the par- 
ticle number density at cell and y,j (= 2;wsin0yAr,A0j) the 
cell volume. The term fi j is the fraction of cell (/, j) which is ion- 
ized, and is set to unity for ionized cells or zero for neutral cells. In 
"boundary" cells it has an intermediate value (see Equation ll9l be- 
low). The particle number density is obtained from the mass density 
by assuming that the gas is entirely neutral hydrogen: «,j = Pij/inu 
(i.e. yu = 1). Essentially the first term on the right-hand side is the 
number of ionizations per unit time in grid column j, and the sec- 
ond term is the number of recombinations per unit time. Thus the 
change in the number of ionized atoms in a single timestep of length 
At is 



Expansion of Ionization front into uniform cloud 



ANj = At 



■sm9jAej'^>-Y,af,jnlV,,j 



(16) 



Having evaluated this number we evaluate the updated number of 
ionized atoms along the column to be 



iVJ+i = N"+ANj. 



(17) 



We then evaluate the number of atoms enclosed at a given radius r, 
along the jlh column as 




Time / yr 



Figure 3. Results of the test simulation, for the expansion of an H ii region 
into a uniform density cloud. Here the number density no = 10''cm"' and 
the ionizing flux <!) = lO^'^s"'. The analytic (Spitzer) solution is plotted 
as a solid black line, with the numerical solutions using two different grid 
resolutions shown as dashed and dotted lines. 



N = V n V 



(18) 



where is the inner grid radius. If A'J^' > N^„c,i then the cell is 
flagged as "ionized", and if N"'^' < A'enc.i then the cell is flagged as 
"neutral". The cell along each ray path where equality is reached is 
flagged as the "boundary cell". In this cell we evaluate the fraction. 



fij, of the cell which must be ionized so that A'"^ 



f.. 



(19) 



In addition, a check is performed to ensure that the recombination 
time-scale at each grid cell which is "ionized", t^^^ = (any', is 
much longer than the dynamic timestep. As long as this condition 
is not violated we can safely neglect radiative cooling of the ionized 
gas. The density of the ionized material is sufficiently small that this 
is not a concern in any of our simulations. 

Following this we have added another subroutine, which ad- 
justs the energy of the gas according to the values the ionization 
flags calculated as above. "Ionized" gas is treated as isothermal 
with a sound speed of lOkms"', "neutral" gas is left unchanged. 
Thus at each cell the energy is adjusted as follows: 



where Chot = lOkms"' 
is evaluated as 



1) 

-1) 



if "neutral" 
if "boundary" 
if "ionized" 



(20) 



, and the sound speed in the boundary cell 



(21) 



In this manner the energy in the boundary cells is evaluated is as the 
weighted mean of the "cold" and "hot" energies. Ccom is evaluated 
as c^^y = Pip on the timestep when a cell first becomes partially 
"ionized". The values of and Cbc are carried to the following 
timestep so that this value can be recovered from Eguation llll if the 
ionization front remains in the same cell for one or more timesteps. 

Lastly, in order to increase the performance of the code we 
have parallelised it, for a shared-memory architecture, using the 



OpenMP formalism''. Tests showed that the numerical results were 
identical to those obtained from running the code on a single pro- 
cessor, and that the parallelisation is around 85-90% efficient when 
running on 4 processors. 



3.2 Testing the code 

To test the accuracy of this algorithm we consider the well-studied 
case of an H ii region expanding into a uniform gas cloud of number 
density no. This problem has a well-known 1-D analytic solution, 
first studied by Soitzei 1 197^. If an ionization front expands spher- 
ically into a uniform gas cloud of density no, then the radius of the 
ionization front n(/) at time t is 

n(0 = rsfl + T— r (22) 



where is the Stromgren radius 



\ Anagnl j 



(23) 



Here ag is the Case B recombination coefficient for atomic hy- 
drogen at lO'^K, which has a value of ag = 2.6 x lO^'^cm's"' 
(Cox 2000). We use a 200 x 200 cell grid covering the range 
r = [OAU, 50AU] and 6 = [0, ?r/2]. Our modified version of zeus2d 
solves this problem numerically, and a comparison between the an- 
alytic and numerical solutions is shown in Fig|3| The code repro- 
duces both the initial Stromgren radius and the power-law rise well, 
accurate at the ^5% level. Runs at higher resolution do not change 
the results significantly, indicating that the procedure is numerically 
converged. We note however that this is not a particularly stringent 
test of the code for the conditions we wish to model, as ionization 
of a disc means that the density profile rises exponentially along 
ray paths. Consequently we expect the ionization front to remain in 
a single grid cell for many more timesteps in a disc simulation than 
in this test, and the treatment of the boundary cell is the primary 

See |http : //www . openmp ■ org| . 
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uncertainty in the computational algorithm. However in a spherical 
geometry the expansion of an H ii region into a cloud with such a 
steeply rising density profile would naturally be Rayleigh-Taylor 
unstable, so this problem cannot be used to test the code. In short, 
we find that our code performs well when considering the expan- 
sion of an Hn region into a uniform density cloud, but note also 
that this is not an especially good test for the problem at hand. 



3.3 Disc model and initial conditions 

In order to model direct photoevaporation it is necessary to set up 
a stable disc configuration in zeus2d. We neglect magnetic fields, 
disc self-gravity and radiation hydrodynamics (these options are 
turned off in the code) and adopt an (r, 8) grid covering a polar an- 
gle of 9 = [0, ;r/2] (i.e. assuming symmetry and simulating only 
one quadrant of the disc). The rotation option in the code, which 
introduces a centrifugal "pseudo-force" is turned on, and acceler- 
ations due to gravity are evaluated using only a point mass placed 
at the origin. The ionization subroutines are included, but initially 
the ionizing flux is set to zero so that they have no effect. Both the 
upper and lower angular boundaries are set to be reflective in order 
to account for the symmetry of the problem, and both the inner and 
outer radial boundaries are set as "outflow" boundaries. We address 
the influence of these boundary conditions later (see Section B31 . 

To determine initia l conditions we use a 1-D reference model 
of the form described bv lClarke etaljj lOOl l . see also Paper II). We 



use the updated wind profile of lFont et al.' r2004', kindly provided 
in numerical form by Ian M cCarthy) rather than the original profile 
of lHollenbach et alJ fl994) used bv lciarke et aljj200ll) . We expect 
direct photoevaporation to become significant once the inner disc 
is drained, so we run the 1-D evolution model to the point where 
the inner disc (inside the draining radius) becomes optically think 
to ionizing photons. The reference model has an ionizing flux of 
<5> = 10'"s"', an initial disc mass of O.OSMq, and an initial accretion 
rate 5.0 x 10^'Moyr"'. At the point when the inner disc is drained 
(at an age of ll.lMyr) the disc has a total mass of O.OOIMq, and 
in order to remove numerical fluctuations we use a functional fit to 
the density profile at this point. The functional fit to the midplane 
density takes the form 

1 



(24) 



This is a R^^ decline at large radii, with a cutoff at radii close to 
some inner edge radius Ri„. p(R,z = 0) oc R^^ is consistent with a 
X oc surface density and constant H/R. This profile peaks at a 
radius Rq, with the parameters related by 



R„ 



(25) 



2(Ro - /?in) 

Thus /?in is the inner truncation radius and Rq the radius at which the 
midplane density peaks. Fitting this profile to the reference model 
gives R,„ = 2.25AU and Rq = 8.25AU, and therefore n = 0.1875. 

Consequently the input parameters to the disc model are as 
follows: M, (stellar mass), Xq (surface density at R = Rg), R,„, Rg, 
H/R and ionizing flux (t>. The initial disc structure is obtained from 
these parameters as follows. The reference midplane density p(R = 
i?o, z = 0) is evaluated as 

So 



p(R = Ro,z = 0) ■■ 



(26) 



where the scale height at Rq is Hq = Rg (H/R). This fixes the con- 
stant of proportionality demanded in Equation 1241 and the cutoff 



index n is evaluated as specified in EQuation l25l Consequently, at 
every point on the (r, 9) grid the mass density and energy density 
are set thus^ : 



p(r, 9) = p(R, z = 0) exp 



2H^I 



(27) 



and 



e(r, 9) 



1 GM, /H\- I z^ \ 

— P(^,. = 0)— (-)exp(-^). (28) 



The functional form is undefined for R ^ Ri„, and so the density 
and energy are set to a constant, small value in this region. (We 
use a value that is 10"'^ times the maximum value.) The radial and 
polar velocities are set to zero. The rotational velocity is set to the 
Keplerian value, with small corrections made to balance the radial 
pressure gradient arising from the choice of density profile. 

Throughout our simulations we adopt a grid which is linearly- 
spaced in both r and 6. There are computational advantages to using 
logarithmic spacing in r (see e.g. . Bate et al..2002.) . in particular the 
fact that the grid cells can all be approximately square (i.e. Ar ^ 
rAd). However a logarithmic grid naturally concentrates the highest 
resolution close to the inner boundary, while we seek to place the 
inner disc edge at larger radius. Also, our simulations do not cover 
a large dynamic range in radius (typically only around one order of 
magnitude). Consequently we find that a logarithmic grid requires 
significantly more grid cells than a linear one in order to achieve 
the same resolution at the inner disc edge, resulting in significantly 
larger CPU requirements, so we adopt a linear grid throughout. 



3.4 Simulations 

As seen above, all of our hydrodynamic disc models are specified 
by the five input parameters M,, So, Ri„, Rq and H/R, plus the ion- 
izing flux O. We adopt M, = IMq throughout. A number of simu- 
lations were run, using the parameters specified in TableQ The ma- 
jority of these simulations use values of Eq , R^ and Rq from the ref- 
erence model. However these parameters were also varied in spe- 
cific cases. Our fiducial model adopts H/R = 0.05 and O = 10*' s"' , 
and we use a suite of simulations to study the elfects of varying 
these parameters. All of the simulations were run for many outer 
orbital times beyond the point where the initial transients died out. 
Here we summarize the simulations conducted: 

Reference This simulation was run as a reference model to 
demonstrate stability of the initial conditions, with the ionizing flux 
= 0. The values of So, ^in and Rq were taken from the reference 
model. We adopt the fiducial value of H/R = 0.05. The grid reso- 
lution was chosen so that Ar = 0.02AU and A9 = With these 
choices the grid cells are approximately square (i.e. Ar ^ rA9) 
at r = 2.5AU. This ensured that the scale-height H was resolved 
into > 5 cells throughout. In order to prevent the run-time becom- 
ing unreasonably long it was necessary to limit the radial range to 
[1.0AU,9.0AU]. 

Edge These simulations were run to study the evolution of the 
inner disc edge. The fiducial model uses the same parameters as 
the reference model above, but with <I> = IC^s"'. Further models 
were run with H/R = 0.1 and O = lO'^^s"' to investigate the effect 
of varying these parameters. 



' Note that the polar and cylindrical coordinates are related by i? = r sin t 
and z = r cos 6. 
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Simulation 




Ns 


\T mmi^ max] 


Sill 




So 


H/R 


<I> 








AU 


AU 


AU 


10"2g cm-^ 




lO't's-' 


Reference 


400 


200 


[1.0,9.0] 


2.25 


8.25 


4.46 


0.05 


0.0 


EdgeI 


400 


200 


[1.0,9.0] 


2.25 


8.25 


4.46 


0.05 


1.0 


Edge2 


400 


200 


[1.0,9.0] 


2.25 


8.25 


4.46 


0.1 


1.0 


Edge3 


400 


200 


[1.0,9.0] 


2.25 


8.25 


4.46 


0.05 


10.0 


ConvTest 


800 


400 


[1.0,9.0] 


2.25 


8.25 


4.46 


0.05 


1.0 


Boundary 


800 


200 


[1.0,17.0] 


2.25 


8.25 


4.46 


0.05 


1.0 


Profile 1 


1200 


100 


[1.0,49.0] 


2.25 


8.25 


4.46 


0.05 


1.0 


Profile2 


1200 


100 


[1.0,49.0] 


2.25 


8.25 


4.46 


0.1 


1.0 


Profiles 


1200 


100 


[1.0,49.0] 


2.25 


8.25 


4.46 


0.05 


10.0 


Profile4 


1200 


100 


[1.0,49.0] 


2.25 


8.25 


4.46 


0.075 


1.0 


LaunchVel 


1200 


100 


[1.0,49.0] 


13.5 


20.0 


10.80 


0.05 


1.0 



Table 1. List of simulations run, showing resolution and physical parameters for each simulation. N,- and Ng are the number of grid cells used in the radial and 
angular dimensions respectively. 



ConvTest This model was run at double the resolution of the 
fiducial model in order to check that the code was numerically con- 
verged. 

Boundary This model was run at the same resolution as the fidu- 
cial model but with double the radial range, in order to study the 
influence of the outer boundary on the results. 

Profile These simulations were run at slightly lower spatial reso- 
lution over a much larger radial range ([1.0AU,49AU]), in order to 
study the mass-loss profile far from the inner disc edge. As with the 
inner edge simulations, 3 simulations were run: the fiducial model, 
H/R = 0.1, and O = 1 0"*^ s" ' . In practice it was found that the mass- 
loss profile was rather sensitive to the disc thickness, and so a fur- 
ther simulation was run with an intermediate value of H/R = 0.075. 

LaunchVel For this simulation the disc was "moved" to a larger 
radius (Rq = 20.0AU), in order to study the launch velocity at the 
ionization front in the regime where the sound speed of the ionized 
gas was greater than the Keplerian velocity. 



4 RESULTS 

The reference model, where the ionizing flux was set to zero, 
demonstrates that the initial conditions are stable. After several 
outer orbital times the only changes to the structure are a slight ver- 
tical expansion of the upper region (due to the failure of the z <ii R 
approximation), and a slight spreading at the inner cutoff radius 
(due to a small pressure mismatch at the cutoff radius: the func- 
tional form has a discontinuity in dP/clR here). In addition there 
is some outflow (and also some spurious reflection) of material at 
the outer radial boundary, due to the artificial pressure gradient in- 
troduced at the boundary. However we are satisfied that the initial 
configuration is stable, and now seek to study the effects of ionizing 
radiation on the disc. 

At this point it is important to note a peculiarity inherent to 
grid-based codes such as zeus2d. When using a grid to specify hy- 
drodynamic variables there are two distinct methods we can use to 
extract data from the simulations. We can either consider the local 
values of the various hydrodynamic variables by looking at their 
values in individual grid cells, or we can study integrated proper- 
ties across many cells. We wish to study the behaviour of the disc at 
the ionization front, as the flow is determined by the hydrodynamic 
properties at the front. As seen in Section|5| this can be done either 
by looking at the density and launch velocity at the ionization front. 




Figure 4. Snapshot of simulation EdgeI at ; = 98yr. Density is plotted 
as a colour scale, with the grid boundaries denoted by solid lines and the 
ionization front by a dashed line. Velocity vectors are plotted at regular 
intervals, but are omitted when they are either smaller than one-fifth the 
length of the reference vector, or when the density is below the minimum of 
the colour-scale. Note, however, that the vectors represent the velocities at 
individual grid cells, and so can be prone to fluctuations. Note also that the 
data is interpolated from the polai' computational grid on to a rectangular 
coordinate grid in order to create the plot, so the resolution of the image 
does not reflect the resolution of the simulation. 



or by considering the integrated mass-loss rates inside fixed radii. 
In terms of the simulations, the density and velocity at the front are 
"single-cell" quantities, whilst the mass outflow rate is a "many- 
cell" quantity. Both are useful, but in general in our simulations we 
find that the "many-cell" outflow rate is more robust against numer- 
ical fluctuations. We also note that our treatment of the ionization 
front inherently results in the measurement of "single-cell" quan- 
tities near the ionization front being resolution-dependent, even if 
the simulation is converged. This is because in order to measure 
properties at the front it is necessary to look at the first "ionized" 
cell: with a different cell size this is at a different distance from the 
front, and therefore always depends on the resolution. 
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4.1 Flow solution 

All of our models show a similar flow field. A snapshot of the fidu- 
cial EdgeI model is shown in Fig|4| The ionization front is per- 
pendicular to the midplane at some inner edge radius, follows the 
disc surface as radius increases and becomes asymptotic to the line- 
of-sight from the ionizing source at large radii. The inner edge ra- 
dius advances slowly with time, with the consequence that the flow, 
while stable, is never entirely steady. The flow is launched almost 
perpendicular to the ionization front. Consequently at large radii 
the flow is launched near to vertically, while at the inner edge the 
launch velocity is radially inwards. Conservation of angular mo- 
mentum prevents the accretion of material through the inner radial 
boundary. The ionized gas at the inner edge follows a ballistic tra- 
jectory, with the closest point of approach visible in the snapshots 
as a sharp decrease in the density of the ionized gas at approxi- 
mately half the inner edge radius (as we would expect for mate- 
rial launched inward at close to the Keplerian speed). Additionally, 
there is a shock front produced by the interaction between the out- 
flowing and inflowing gas near to the inner disc edge, which again 
is clearly visible in the snapshots as a discontinuity in the density of 
the ionized gas. However the flow between the disc surface and the 
shock front is supersonic, so the shock front does not influence the 
flow at the ionization front. As the flow progresses along stream- 
lines it becomes almost radial, but with streamlines which diverge 
from the inner disc edge rather than the origin. 

In general the sonic surface is very close to, but slightly above, 
the ionization front (see Fig|5}, suggesting that the flow is launched 
sub-sonically before rapidly becoming supersonic, as we expect for 
any steady wind. However we note at this point that the velocity 
field near to the ionization front is not always stable. This is because 
the ionization front is not coincident with the grid columns. Conse- 
quently, as the front evolves across cell boundaries it can produces 
"kinks" in the front, which momentarily result in converging flows 
near to the front and are responsible for the striations visible in the 
flow snapshots. As a result the velocity structure near to the front 
can be rather messy, while remaining smooth far from the front. As 
long as the ionization front remains unresolved (as assumed in the 
algorithm) this problem will persist. In practical terms this means 
that this effect is be present in all our simulations, regardless of the 
resolution. 

A further feature of the algorithm is the treatment of material 
which flows in the 0-direction, across different ray paths. In our al- 
gorithm, the change in the number of ionized atoms is evaluated in 
such a way that any material flowing into a ray path is assumed to be 
neutral (see Eguations llSHlSL Thus flow across the front from the 
neutral to the ionized region is treated correctly. However at higher 
latitudes, in the ionized region, this is incorrect and results in the 
ionization front being placed at too small a radius in some cases, re- 
sulting in some rather odd-looking front geometries at early times. 
However in practice this is a transient efi'ect which dies out as the 
flow becomes more stable, and only affects the simulations at early 
times. Consequently we do not deem it to be a significant problem, 
as a similar time is needed for the other initial transients to disap- 
pear. 

4.2 Convergence, boundary effects and other numerical 
issues 

The first computational consideration is the possible influence of 
numerical effects on the results obtained from the simulations. Two 
specific issues are addressed: numerical convergence, and the influ- 




Figure 5. Snapshots of simulation EdgeI at / = 98yr & 291yr. The plot is as 
in Fig|4] but with the sonic surface plotted instead of velocity vectors. The 
sonic surface is shown as a solid line, with tick marks indicating the sub- 
sonic side of the surface. Note the "messy" velocity structure near to the 
ionization front in the second snapshot. Comparison of the two snapshots 
also highlights the influence of the outer boundary condition on the outer 
region of the unperturbed disc. By ? = 291yr (=: 10.8 orbital periods at the 
outer boundary) neutral gas has begun to "fall off' the disc at large radii due 
to the artificial pressure gradient introduced by the boundary condition, as 
can be seen by comparing the two snapshots. 



ence of the artificial boundaries imposed upon the simulations. In 
order to address the issue of numerical convergence we performed 
the ConvTest simulation. This simulation is identical to the fidu- 
cial model (EdgeI), but has double the spatial resolution in both 
dimensions. In terms of the disc scale-height, HjrbS = 6.3 in sim- 
ulation EdgeI, and HIrM) = 12.7 in ConvTest. To test conver- 
gence we evaluated the mass flux through a series of surfaces at 
fixed r in the ionized region, multiplying by a factor of two to ac- 
count for flow from both sides of the disc. The surfaces used were 
not near to either the inner disc edge or the computational bound- 
ary in order to minimize their effects: surfaces at r = 5, 6, 7, & 
8AU were used. The quantity M(< r) should be independent of 
resolution, and so it is instructive to study the ratio of the different 
mass-fluxes between the two simulations. The ratio of mass-fluxes 
(EdgeI/ConvTest), time-averaged over f = to / = 80yr (approx- 
imately three outer orbital periods), was found to be 1.04 ± 0.06, 
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1.02 ± 0.05, 1.04 ± 0.05 & 1.06 ± 0.06 at radii of 5, 6, 7 & 8AU 
respectively. These results suggest that the solution is converged, as 
the discrepancies between the two simulations are smaller than the 
fluctuations within each individual simulation. The results also in- 
dicate that the numerical scheme is accurate to around ±5%. Runs 
at lower resolution suggest that the solutions remain converged to 
better than 20% accuracy as long as the disc scale-height H(R) is 
revolved into at least 3 grid cells, and that this accuracy improves 
as the radial range increases. (This is due to Ar being fixed: a ra- 
dial length of H is resolved into more grid cells at larger radius) 
Thus for H/R = 0.05 the minimum resolution requirement for 
a marginally converged solution is A9 = n/200 (i.e. 100 cells in 
the angular direction). We adopted double this resolution wherever 
possible, but were forced to use this minimum resolution when con- 
sidering a large radial range in order keep the CPU requirements to 
a reasonable level. Consequently we are satisfied that the simula- 
tions are numerically converged, although we note that the numer- 
ical accuracy of the Profile simulations is likely only ±15%. 

The next issue to address is the influence of the computa- 
tional boundaries on the simulations. Both of the angular bound- 
aries are axes of symmetry, and so the use of reflective boundaries 
here is physically reasonable. Further, zeus2d is exact in its treat- 
ment of reflective boundary conditions, so they are treated correctly 
by the code. We note that the angular momentum barrier is suffi- 
ciently large that, initial transients aside, no mass actually reaches 
the R = angular boundary. Both of the radial boundaries are set 
to be outflow boundaries. zeus2d is exact in its treatment of out- 
flow boundaries as long as the flow is supersonic and along grid 
columns. However sub-sonic outflow is known to produce spuri- 
ous reflection at the boundaries, and supersonic outflows which are 
not along grid columns can also be problematic. We note than in 
our simulations accretion through the inner boundary is negligible 
(and is exactly zero aside from initial transients), and so the inner 
boundary condition has no effect on the results. 

In order to study the influence of the outer radial boundary 
condition we performed the simulation Boundary. This simula- 
tion has exactly the same parameters and resolution as the fiducial 
model (EdgeI), but covers double the radial range. Consequently it 
is possible to study the influence of the outer boundary in the fidu- 
cial model by comparing the region near to the boundary in simula- 
tion EdgeI to the same region in simulation Boundary (where it is 
in the middle of the computational domain). We evaluated the mass 
outflow rates of both ionized and cold material through surfaces 
at r = 7, 8, & 9AU in both simulations. The effect of the bound- 
ary on the ionized material is to decelerate the flow somewhat near 
to the boundary: the flow rates at 7AU and 8AU are essentially 
equal in the two simulations (time-averaged ratios of 1.02 ± 0.06 
and 1.01 ± 0.06), but the outflow rate is somewhat smaller at 9AU 
in EdgeI than in Boundary (0.89 ± 0.06). We hypothesise that this 
is caused by spurious reflections at the boundary, due to the flow di- 
rection not being radial (as supersonic outflow along grid columns 
should be treated exactly). However we overcome this issue simply 
by neglecting the region near to the boundary: in the following sec- 
tions we do not evaluate any variables at cells which are in the outer 
10% of the grid. The efl"ect of the boundary on the neutral material 
is to introduce a radial pressure gradient across the outer boundary, 
leading to some outflow of material. The region influenced by this 
grows with time (as seen in Fig|5}, and is also larger in simulations 
with greater H/R. This effect is not significant in most of the sim- 
ulations, but care must be taken when analyzing simulations with 
large H/R or which have run to very large problem time. 



Simulation 


Time Interval 


Best-fitting 


parameters 




yr 


Power-law index b 


Normalisation C 


EdgeI 


15-1000 


-1.54 


0.56 


Edge2 


15-500 


-1.48 


0.79 


Edge3 


10-500 


-1.46 


0.64 


Mean 




-1.49 


0.66 



Table 2. Best-fitting parameters to inner edge density. The parameters are 
fit using a simple least-squares algorithm. The final row shows the mean 
values from the three simulations. The typical (Icr) uncertainties in the fits 
are ±0.1 in the power-law index and ±25% in the normalisation constant. 



4.3 Inner edge density 

Three simulations were run in order to study the evolution of the 
inner edge of the disc. Our simple analytic argument (Equation^ 
suggests that the density at the ionization front at the disc midplane 
should scale as 

n,n = C\ ■-] , (29) 

\4na(H/Rk„Rl) 

where C is an order-of-unity constant. The scaling constant ac- 
counts for the fact that the radial scale-height is not exactly equal 
to the vertical scale-height, and also for the attenuation of the radi- 
ation field at R < Ri„. We now seek to compare this prediction to 
the results of the simulations. 

A problem arises when considering the density at the ioniza- 
tion front in the simulations. The ionization front is not resolved, 
and so the density must be measured in the grid cell adjacent to the 
boundary cell. This value decreases as the front advances across 
each cell, and so for the same value of Ri„ (i.e. location of the 
boundary cell) a spread of densities is observed. However by av- 
eraging over a large number of timesteps the noise introduced is 
minimized. Additionally, as mentioned in Section^above, the eval- 
uation of "single-cell" quantities is always resolution dependent. 
Therefore, while it is possible to fit the power-law index accurately 
we find that the value of the normalisation constant C depends 
somewhat on the grid resolution. 

In order to test Equation above we performed a fit to the 
data, using the form 

O 



Pin = Coth 



1/2 



Rt 



(30) 



[4na(H/R)j 

A two-parameter least-squares fit was performed to find the best- 
fitting power-law index b and normalisation constant C. The results 
of these fits are shown in Table|5|and Fig|6| For simulations Edge2 
and Edge3 a smaller range in problem time was used in order to 
lessen the influence of the outer boundary. (The outer boundary is 
more significant in both of these simulations that in EdgeI, due to 
the larger H/R in Edge2 and the faster progression of the inner edge 
in Edge3.) Due to the relatively small change in R,„ over the course 
of the simulations there is a degeneracy between the two param- 
eters, leading to rather large error bounds. The best-fitting values, 
averaged over all three simulations, are b = -1.49 and C = 0.66. 
The typical (Icr) uncertainties in the fits are ±0.1 in b and ±25% in 
C. Thus the results are entirely consistent with a power-law index 
of -1.5, as predicted analytically. The ConvTest simulation does 
not cover a sufficient radial range in R^y, to constrain the power-law 
index well, giving a best-fitting value of b = 1.7 ± 0.4. However 
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Inner edge density: *=10*'s \ H/R=0.05 




2 3 4 5 6 



Ri„ / AU 

Figure 6. Inner edge density plotted as a function of edge radius for simula- 
tion EdgeI. Points are plotted eveiy lyr of problem time, from 15-lOOOyr 
A power-law fit is shown as a dashed-line, with index b = - \.5 and normal- 
isation constant C = 0.6. 

when the index is set lob = -1.5 the best-fitting value of the nor- 
malisation constant is 0.83. Thus it seems that the value of C is in- 
deed somewhat resolution-dependent, but seems to lie in the range 
0.5-1.0. Note also that the determination of C to be of order unity 
verifies our earlier assumption (Section|5J that the thickness of the 
ionized layer is comparable to H, at least along the disc midplane. 

4.4 Launch velocity 

A number of simulations were run to try to constrain the launch 
velocity and determine the value of the constant D in Equation |5| 
However, as mentioned above, the velocity field near to the ion- 
ization front fluctuates and it was not possible to obtain good con- 
straints on the launch speed. In general it seems that the launch 
speed is approximately constant, as there is no significant variation 
with radius. However the fluctuations in the velocity at any given 
radius are large. In general the launch speed seems to be slightly 
sub-sonic, although it can become supersonic at times. The best- 
fitting values of D show a large scatter, and cannot be constrained 
better than to say that the launch speed is, as expected, of order the 
sound speed. It is necessary to consider "many-cell" variables in 
order to obtain a more accurate result. 

4.5 Wind profile 

In order to study the wind profiles produced by the model we con- 
ducted the Profile series of simulations. These simulations cover 
a large radial range, enabling the study of the wind profile over a 
wider range of radii than possible with the higher resolution edge 
simulations. The method adopted to analyse the results of these 
simulations was to study the mass outflow rate across given radial 
shells at each timestep. At radii outside the inner disc edge the out- 
flowing ionized mass through a shell of constant r comes from the 
disc surface inside a cylindrical radius R = rsin Of (where Of is the 
angular coordinate of the ionization front). Consequently by sam- 
pling M(< r) at many radii it is possible to study the mass-loss 
profile t(R). 

In order to study the simulations in more detail we compared 
our analytic prediction for M(< R) (derived in Section |2j with 



the values obtained from the simulation. In general we find that 
a power-law form for the shape function 

fix) = x-" (31) 

provides a good fit to the profile, although the values of a varies for 
different values of HIR (i.e. we fit a different shape function for dif- 
ferent H/R). Substituting this form into the the analytic expression 
for M(< R) fEQuation ll3> and integrating gives 




for a i= 2. Therefore by comparing to the numerical results it is pos- 
sible to obtain values for the combined scaling constant (CD) and 
the power-law index a. Our procedure for fitting these parameters 
is as follows. 

At each timestep we evaluate M(< R) at a large number of 
different radii on the grid. The mass outflow rate is evaluated as 

M(< r;) = 2 ^ Pijfrd.j) {2nAij sin 6*^ A6)) (33) 

J 

in the ionized region. Here Kr(i,j) is the radial component of the fluid 
velocity at cell (/, j), the last term represents the area of the outer 
radial surface of each cell, and the initial factor of 2 accounts for 
flow from both sides of the disc. The smallest value of r, used was 
0.5AU larger than the inner edge radius at that timestep, and we 
do not evaluate M(< Vj) in the region where the outer boundary 
becomes significant (the outer 10% of the grid for most models, 
but somewhat larger where HIR > 0.05). We then solve for the 
values of a and {CD) which give the best-fit to a straight line in the 
M(< /?)/[l -{RmIRT ^\ plane. These values fluctuate somewhat 
with time, but by studying a large number of timesteps it is possible 
to obtain best-fitting values. In each case it is necessary to wait 
until the simulation stabilises before studying the mass-loss profile: 
typically this takes around one outer orbital time. The evolution of 
the parameters with time for the Profile 1 simulation are shown in 
FigEl and ths mean values of the constants for each simulation are 
shown in Table|3| 

In general we find that a single power-law index and a sin- 
gle normalisation constant provide a good fit to the mass-loss pro- 
file for fixed HIR. For HIR = 0.05 the best-fitting values are 
a = 2.42 ± 0.09 and (CD) = 0.235 ± 0.02: examples of this fit 
are shown in Fig|8| The single value obtained confirms that the 
mass-loss rate scales as <1>''^, which in turn confirms th at the flow 
is "recombination-limited" (see lHoUenbach et"al .11994 . Note also 
that this value of (CD) is consistent with C - 0.6, as estimated 
from the inner edge density, if D ^ 0.4. This suggests that the 
typical launch velocity is around 0.4 times the sound speed, a value 
consistent with numerical simulatio ns of the photoeva porative wind 
driven by the diffuse radiation field front et alj2'004l) . 

However the power-law index is rather sensitive to the value 
of HIR, as seen in Tabled This is not entirely surprising, as the 
disc density decreases exponentially with (zlH)^, meaning that the 
density at which ionization balance occurs is rather sensitive to the 
disc thickness. This obviously has a knock-on effect for the inte- 
grated mass-loss rates, with thicker discs resulting is less eflicient 
mass-loss at radii beyond the inner disc edge. However while the 
shape of the mass-loss profile in terms of t.(R) varies significantly 
with the value of HIR, the elfect on the integrated mass-loss rate 
is not especially strong. This is illustrated in Fig|9| which shows 



Best-fitting power— law Index for simulation PR0FILE1 



Disc photoevaporation I: hydrodynamic models 

Best— fitting normalisation constant for simulation PR0FILE1 



11 




500 600 700 800 900 1000 500 600 700 800 900 1000 

t / yr t / yr 

Figure 7. Best- fitting parameters to tlie wind profile as a function of time for simulation Profile 1. The left-hand panel shows the power-law index a, the 
right-hand panel the normalisation constant (CD). In both plots the mean value is shown as a dashed line, with the ±1(T error bounds shown as dotted lines. 



Simulation 


Time Interval 


H/R 


D 


Time-averaged parameters 




yr 




lO^'s-l 


a 


(CD) 


Profile 1 


500-1000 


0.05 


1.0 


2.33 ± 0.09 


0.21 ±0.02 


Profile2 


250-500 


0.1 


1.0 


4.50 ±0.13 


0.60 ± 0.02 


Profiles 


250-500 


0.05 


10.0 


2.49 ± 0.08 


0.26 ± 0.02 


Profile4 


250-500 


0.075 


1.0 


3.38 ±0.18 


0.43 ± 0.03 



Table 3. Time averages of the scaling parameters for the mass-loss profile, evaluated from each simulation, a is the power-law index, and (CD) is the 
normalisation constant. The errors quoted are (Icr) standard deviations. 



the best-fitting mass-loss profiles for three different values of H/R. 
Larger values of H/R result in mass-loss profiles which are rather 
more concentrated towards the inner disc edge, but the total inte- 
grated mass-loss rates differ only by a factor of 2-3. 



5 DISCUSSION 

In order to construct these models we have made a number of sim- 
plifications and approximations, and now seek to address their sig- 
nificance. The most important of these is the use of the on-the-spot 
(GTS) approximation to simplify the radiative transfer problem. 
Throughout the simulations we have adopted Case B recomb ina- 
tion coefficient, which has a value of = 2.6 x 10^'^^cm^^s 
12000 !). Adopting the Case A value (^a = 4.2 x lO^'^cm^s"') has a 
negligible effect on the results, as the density in the ionization bal- 
ance equation depends on a"''^. However we must also consider 
the possibility that recombination photons can be absorbed non- 
locally: were this to be significant the GTS approximation would 
fail. In the simulations the ionized region is always optically thin to 
ionizing photons. (In fact, the density in the ionized region is suffi- 
ciently low that it would not be optically thick to ionizing photons 
even if the gas were entirely neutral.) Similarly, the density in the 
cold disc is sufficiently large that the mean free path of an ionizing 
photon is always much smaller than the length of a grid cell. Conse- 
quently, if we consider recombinations at a given point on the front 
we see that any photons emitted in an upward direction escape the 
system, whereas any photons emitted downwards are absorbed lo- 



ntegrated mass-loss profiles: *=10 s , R| =3.5AU 




R / AU 

Figure 9. Best-fitting mass-loss profiles for various values of H/R. The 
profiles are evaluated for fixed values of <I) = 10'*'s"' and Rj,, = 3.5AU, 
using the best-fitting parameters listed in Table l3l The solid line shows the 
profile for H/R = 0.05, the dotted line H/R = 0.075, and the dashed line 
H/R = 0.10. While the form of the profiles varies somewhat with H/R, the 
integrated mass-loss rates are not especially sensitive to the disc thickness. 



cally. If we look at the geometry of the ionization front (see Figs.|4| 
&|5j, we see that the front approaches the line-of-sight from the 
source asymptotically at large radii. Therefore there is only a very 
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ntegrated mass-loss profile: *=10*'s \ H/R=0.05, t=750yr Integrated mass-loss profile: *=10''^s \ H/R=0.05, t=350yr 




10 20 30 40 10 20 30 40 



R / AU R / AU 

Figure 8. Snapshots of the mass-loss profile from the simulations. The left-hand panel shows the results for simulation ProfileI (<I) = lO'^'s"', HjR = 0.05) 
att = 750yr, the right-hand panel simulation Profiles (d) = lO^^s"', H/R = 0.05) att = 350yr. In both plots the dots are the data from the simulations, with 
the solid line showing the analytic form of the mass-loss profile (Eauation ll3l for a = 2.42 and (CD) = 0.235: in both cases the fit is extremely good. 



small solid angle which peimits recombination photons to be ab- 
sorbed elsewhere along the front, and their influence is negligible 
by comparison to the direct field. Consequently we are satisfied that 
the OTS approximation is valid for this geometry. 

A further simplification related to the the use of the OTS ap- 
proximation is that recombination photons produced in the flow 
region are similarly neglected. This diffuse field is th e field that 
drives the disc wind in the models of Hollenbach et ^ and 
so by comparison to their work it is possible to estimate its con- 
tribution. Hollenbach et al. ( 1994) found that the diffuse field pro- 
duces around 10% of the flux of the direct field in the static region. 
Here the diffuse field must come from ionized material that is flow- 
ing away from the disc, and so may be somewhat less efficient than 
in the static case. However if we neglect the motion of the flow it 
seems that neglecting the diffuse field produced by recombinations 
in the flow results in our model under-estimating the radiation field 
at the front by around 10%. We suggest, however, that a more de- 
tailed radiative transfer calculation is needed in order to confirm 
this. We note that both of the simplifications of the OTS approx- 
imation result in our model under-estimating the ionizing flux at 
the ionization front, and therefore suggest that our mass-loss rates 
can be considered as lower-limits to the true mass-loss rates due 
to direct photoevaporation. Moreover the wind rate scales as O''-, 
and so a small increase in the effective value of O does not have a 
strong influence on the total mass-loss rate. 

The second simplification we have made is the equation of 
state adopted in the "boundary" cells (Equations 1201 & This 
approximation may be significant, as the treatment of the boundary 
cell can have a significant effect on the flow solution. As mentioned 
above, we treat the boundary cells in such a way that any material 
flowing "laterally" is assumed to be neutral. In practice the pressure 
gradients across the ionization front are upwards, from "neutral" to 
"ionized", and so this should predict the correct behaviour. How- 
ever with no model with which to compare our results it is difficult 
to quantify the accuracy of this approximation. As seen in Sec- 
tion comparison to the Spitzer solution suggests an accuracy 
of around 5%, but we also note this is not an especially good test 
of the code's applicability to the direct photoevaporation problem. 



Comparison of our results with independently obtained solution in 
the future will be a valuable test of the accuracy of the code. 

As noted in Section |2| our model neglects any possible dust 
opacity between the star and the inner disc edge. The presence of 
dust in the "inner hole" could attenuate the flux of ionizing pho- 
tons reaching the disc, and consequently reduce the wind rate pro- 
duced by direct photoevaporation. However it is not clear whether 
or not the viscous accretion of the inner gas disc will also remove 
the dust, and a detailed treatment of the gas-grain dynamics during 
the transition phase is beyond the scope of this work. S ome in- 
sight can be gained from the work of'Take uchi et alj i2005h . which 
considers the differential motions of dust and gas in a simple pho- 
toevaporating disc model. In general they find that dust is rapidly 
accreted with the gas when the inner disc drains, but there is a com- 
pl icated dependence on the various model parameters. The results 
of ^Takeuc hi et alj j20o3) therefore suggest that the effect of dust 
opacity in the inner part of the disc will be small, but we note that 
this issue warrants further study. 

As seen above, the mass-loss profiles obtained by the model 
are rather sensitive to the disc scale-height. We note also that 
real TT discs are expected to show significant "flaring" (e.g. 
Kenvon & Hartmann 1987), and are not generally consistent with a 
constant H/R ratio. A flaring disc would show a slower decrease in 
the incidence term sin/3 (Equation|5j with R than a disc with con- 
stant H/R. Consequently a flaring disc is expected to show greater 
mass-loss at large radii than seen in our models. However once 
again this suggests that our models provide a robust lower-limit 
to the mass-loss rates produced by direct photoevaporation. 

Lastly, it is instructive to compare the mass-loss profile we 
obtain to that predicted by diffuse photoevaporation. In the diffuse 
model the wind profile takes the form 

Swind R-'/' , (34) 

with the profile normalised at R„ regardless of the instantaneous 
location of the inner disc edge (Hollenbach et al. 1994). Therefore 
as the inner disc edge moves outward to Ri^ > Rg, the integrated 
mass-loss rate decreases significantly. However when direct pho- 
toevaporation is considered we find that while the decline in the 
radial power-law is similar (for H/R = 0.05), the profile is instead 
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normalised at the instantaneous inner disc edge (see Eguation llOt . 
As discussed in Section |2| th e form of the direct wi nd is qualita- 
tivehy similar to that found by HoUenbach et for the case 

of a strong stellar wind. This form results in a significantly larger M 
as the inner disc edge evolves outwards, as the integrated mass-loss 
rate scales wi\hR\l^ (Eguation lOl . Additionally, the more effective 
radiative transfer in the case of direct photoevaporation results in a 
mass-loss rate that is larger than that expected in the diffuse case, 
by a factor of ^ 5-10 at ~ R^. Therefore we expect that direct 
photoevaporation is much more efficient than diffuse photoevapo- 
ration in dispersing the outer part of the disc. This provides a pos- 
sible sol ution to the "oute r disc problem" of the original UV-switch 
model I Clarke et al ]!2001). We explore the effect of this modified 
mass-loss profile on the evolution of TT discs in a companion paper 
(Paper 11). 



6 SUMMARY 

We have highlighted a significant flaw in disc evolution models 
which incorporate photoevaporation, namely that such models ne- 
glect the direct radiation field after the inner disc has been drained. 
We have modelled the photoevaporative wind produced by the di- 
rect field, firstly using a simple analytic treatment and then using 
detailed numerical hydrodynamics. Our models show that once the 
inner disc has drained the wind due to the direct field dominates 
over that due to the diffuse field. The presence of such a wind leads 
to a significantly faster dispersal of the out er disc tha n seen in ex- 
isting models of disc photoevaporation (e.g. lClarke et al..200D . as 
the integrated mass- loss rate scales with Rlp as i?i„ grows (cf. ^j^''^ 
in the diffuse case, iHollenbach etalJll994 . We suggest that this 
may be a solution to the "outer disc problem" in these models. Our 
analysis results in a simple form for the mass-loss due to direct pho- 
toevaporation, and we explore the effects of this wind on models of 
disc evolution in a companion paper (Paper II). 
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